Preparations

Load the necessary libraries

library(mgcv)      #for GAMs
library(gratia)    #for GAM plots
library(broom)     #for tidy output#
library(emmeans)   #for marginal means etc
library(MuMIn)     #for model selection and AICc
library(tidyverse) #for data wrangling
library(DHARMa)    #for simulated residuals
library(performance) #for residual disagnostics
library(see)        # to visualize residual diagnostics

Scenario

In a chapter on time series analysis, Reed et al. (2007) presented Hawaiian longitudinal waterbird survey data. These data comprise winter counts of various species of stilts, coots and moorehen along with year and the previous seasons rainfall. Here, we will explore the temporal patterns in the Kauai Moorhen.

Moorhen

Format of reed.csv data file

Year Stilt.Oahu Stilt.Maui Coot.Oahu Coot.Maui Moorhen.Kauai Rainfall
1956 163 169 528 177 2 15.16
1957 272 190 338 273 NA 15.48
1958 549 159 449 256 2 16.26
1959 533 211 822 170 10 21.25
1960 NA 232 NA 188 4 10.94
1961 134 155 717 149 10 1 9.93
Year - a continuous predictor
Stilt.Oahu - the abundance of the Oahu stilt
Stilt.Maui - the abundance of the Maui stilt
Coot.Oahu - the abundance of the Oahu coot
Coot.Maui - the abundance of the Maui coot
Moorhen.Kauai - the abundance of the Kauai moorhen
Rainfal - the number of centimeters (or inches) of rain

Read in the data

reed = read_csv('../public/data/reed.csv', trim_ws=TRUE)
## Rows: 48 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (7): Year, Stilt.Oahu, Stilt.Maui, Coot.Oahu, Coot.Maui, Moorhen.Kauai, ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(reed)
## Rows: 48
## Columns: 7
## $ Year          <dbl> 1956, 1957, 1958, 1959, 1960, 1961, 1962, 1963, 1964, 19…
## $ Stilt.Oahu    <dbl> 163, 272, 549, 533, NA, 134, 175, 356, 485, 184, 242, 20…
## $ Stilt.Maui    <dbl> 169, 190, 159, 211, 232, 155, 282, 170, 164, 162, 253, 1…
## $ Coot.Oahu     <dbl> 528, 338, 449, 822, NA, 717, 12, 169, 98, 112, 77, 106, …
## $ Coot.Maui     <dbl> 177, 273, 256, 170, 188, 149, 205, 108, 79, 53, 75, 80, …
## $ Moorhen.Kauai <dbl> 2, NA, 2, 10, 4, 10, 12, 10, 8, NA, 17, 7, 44, 50, 26, 1…
## $ Rainfall      <dbl> 15.16, 15.48, 16.26, 21.25, 10.94, 19.93, 12.63, 20.09, …

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{Pois}(\lambda_i)\\ log(\lambda_i) =\beta_0 + f(Year_i) + f(Rainfall_i) \]

where \(\beta_0\) is the y-intercept. \(f(Year)\) and \(f(Rainfall)\) indicate the additive smoothing functions of Year and Rainfall respectively.

Refit and validate models

Partial plots

draw

draw(reed.gam3)

draw(reed.gam3, residuals = TRUE, scales = 'free') 

plot

#plot(reed.gam3, pages=1)
#plot(reed.gam3, pages=1, shift=coef(reed.gam3)[1])
#plot(reed.gam3, pages=1, shift=coef(reed.gam3)[1], scale=0)
plot(reed.gam3, pages=1, shift=coef(reed.gam3)[1], trans=exp,
     resid=TRUE, cex=4, scale=0)

Model investigation / hypothesis testing

summary

summary(reed.gam3)
## 
## Family: Negative Binomial(4.078) 
## Link function: log 
## 
## Formula:
## Moorhen.Kauai ~ s(Year, k = 20, bs = "cr") + s(Rainfall, bs = "cr")
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.81247    0.08035   47.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df  Chi.sq p-value    
## s(Year)     7.361  9.070 164.116  <2e-16 ***
## s(Rainfall) 1.452  1.754   0.475   0.759    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.796   Deviance explained =   81%
## -REML = 214.57  Scale est. = 1         n = 45

Conclusions:

  • the average number of Moorhen is 3.81 on the link scale or 45.26 on the response scale. This number corresponds to the average number of Moorhens expected for the average year with the average rainfall.
  • there is evidence of a signficantly wiggly change in Moorhen numbers over time.
  • there is no evidence of a wiggly rainfall trend
  • we might consider dropping the smoother for rainfall in preference for a regular linear parametric term.

tidy

tidy(reed.gam3)

Conclusions:

  • the average number of Moorhen is 3.81 on the link scale or 45.26 on the response scale. This number corresponds to the average number of Moorhens expected for the average year with the average rainfall.
  • there is evidence of a signficantly wiggly change in Moorhen numbers over time.
  • there is no evidence of a wiggly rainfall trend
  • we might consider dropping the smoother for rainfall in preference for a regular linear parametric term.

Further modelling

Summary figures

References

Reed, J. M., Elphick C. S., Zuur A. F., and Ieno E. N. 2007. “Analysing Ecological Data.” In, edited by Zuur A. F., Ieno E. N., and Smith G. M., 615–32. Springer, New York.